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We study how the interplay between the memory immune response and pathogen mutation affects 
epidemic dynamics in two related models. The first explicitly models pathogen mutation and in- 
dividual memory immune responses, with contacted individuals becoming infected only if they are 
exposed to strains that are significantly different from other strains in their memory repertoire. The 
second model is a reduction of the first to a system of difference equations. In this case, individuals 
spend a fixed amount of time in a generalized immune class. In both models, we observe four fun- 
damentally different types of behavior, depending on parameters: (1) pathogen extinction due to 
lack of contact between individuals, (2) endemic infection (3) periodic epidemic outbreaks, and (4) 
one or more outbreaks followed by extinction of the epidemic due to extremely low minima in the 
oscillations. We analyze both models to determine the location of each transition. Our main result 
is that pathogens in highly connected populations must mutate rapidly in order to remain viable. 



PACS numbers: 05.45.-a, 



.75.-k, 87.23.Ge, 87.19.Xx 



I. INTRODUCTION 

The memory immune response enables humans and 
other animals to rapidly clear, or even prevent altogether, 
infection by pathogens with which they have previously 
been infected. For example, we typically contract chicken 
pox only once in our lifetime because of the effectiveness 
of the memory immune response, and vaccines are de- 
signed around the knowledge that our immune systems 
will more efficiently fight foreign invaders if already ex- 
posed to something very similar. Consequently, it is easy 
to imagine why some pathogens, such as influenza, use 
a strategy of disguise to survive in a host population. 
In most cases, this disguise is facilitated by mutation: 
pathogens permanently change their genetic content in 
order to alter their appearance to the host immune sys- 
tem. With enough mutations, a pathogen will ultimately 
be unrecognizable to the immune system of a host that 
has previously been infected with one of its ancestors. In 
this paper, we study the dynamics of populations that 
lose immunity via this route. 

In epidemiological models, host populations are tra- 
ditionally categorized into three states: susceptible to 
infection (S 1 ), infected (/), and removed or immune (R). 
The succesion of states we will study is depicted in the 
following diagram: 
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Models that describe such an epidemiological cycle are 
referred to as SIRS models. While there is a vast lit- 



erature covering models in which the "loss of immunity" 
step is not considered (referred to as SIR models; see 
for example the classic texts by Bailey 0] , Anderson and 
May H and the recent review by Hethcote ||), compar- 
atively little work has been done to understand how the 
nature of the R — > S transition affects the dynamics of 
an epidemic. 

In principle, the transition depends on the strain to 
which one is exposed (the challenge strain), in addition 
to one's previous history of infection. We thus begin our 
analysis with a computational "bitstring model" in which 
different pathogen strains are represented by bitstrings 
that can mutate. In this model, immunity depends ex- 
plicitly on the history of strains with which one has been 
infected. We find four fundamentally different types of 
behavior, depending on parameters: (1) pathogen extinc- 
tion due to lack of contact between individuals, (2) en- 
demic infection (steady state infection), (3) periodic epi- 
demic outbreaks (sustained oscillations), and (4) one or 
more outbreaks followed by extinction of the epidemic 
due to extremely low minima in the oscillations ("dy- 
namic extinction"). 

We then develop a difference equation model in which 
the nature of immunity is significantly simplified. Instead 
of acquiring indefinite immunity to a specific pathogen, 
individuals in this reduced model spend a fixed num- 
ber of time steps in a generalized immune class be- 
fore being returned to the susceptible population. This 
SIR1R2 ■ ■ ■ RnS model has been studied extensively by 
Cooke et al. Q and also by Longini || , who used stochas- 
tic methods to investigate the model with immunity last- 
ing only a single time step. The same four qualitatively 
different dynamics seen in the bitstring model are also 
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observed for this model. We extend Cooke's work by de- 
riving the location of the onset of oscillatory dynamics 
in any dimension (which is determined by the number of 
time steps spent in the immune class). Stability is lost 
through a Hopf bifurcation, and changes in the model 
parameters can increase the resulting limit cycle ampli- 
tude to the point that the minimum becomes extremely 
small. We establish a criterion for "dynamic extinction" 
(for which the minimum fraction infected in the limit cy- 
cle oscillation is less than 1 /N, where N is the population 
size) and construct an asymptotic approximation for the 
location of this extinction transition. 



II. BITSTRING MODEL 

The immune response recognizes foreign molecules in 
a highly specific way. Individual immune cells or anti- 
bodies that recognize a protein from one pathogen strain 
may be unable to recognize a similar protein derived from 
another pathogen strain. Because protein sequence and 
structure is determined by the genetic content of an or- 
ganism, immune responses to a pathogen are in fact spe- 
cific to the pathogen's genetic content. We thus introduce 
a bitstring model in which pathogen strains are repre- 
sented by bitstrings, where the bitstring is regarded as 
an abstract representation of a pathogen's genetic code. 
Hosts are immune to infection by pathogens that are 
highly similar to pathogens with which they have been 
previously infected, where similarity between two strains 
is measured by hamming distance H . 

The model consists of N individuals who are either un- 
infected or infected with a strain represented by one of 
the 2 e possible bitstrings of fixed length I . Individuals 
keep a record of all the strains with which they have been 
infected, and we refer to these histories as their memory 
repertoire. Once per time step, each infected individual 
exposes z others by selecting individuals uniformly at 
random from the entire population [gj. The susceptibil- 
ity of a contacted individual is determined by comparing 
the bitstring of the challenge strain with the bitstrings 
of all strains in the memory repertoire. Specifically, an 
individual is susceptible if h m i n > /ithr, where /i t hr is a 
parameter, and h m { n is the smallest hamming distance 
between the challenge strain and any strain in the in- 
dividual's memory repertoire. With probability \x the 
challenge strain mutates by flipping one randomly cho- 
sen bit; otherwise the strain remains unchanged. In each 
case considered here, we keep \x fixed at 0.1 and vary the 
threshold hamming distance /ithr! these two parameters 
are inversely proportional. Infection lasts a single time 
step, and strain transmissibility is determined exclusively 
by individual immune responses, i.e., in an entirely sus- 
ceptible population no strain is more fit than any other. 

In the first time step, one randomly chosen individual 
is infected with the bitstring 000 • • ■ 000 and all others 
have never been infected. From this initial condition we 
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FIG. 1. The dynamics of infection in the bitstring model, 
(a) Fraction of population infected versus time for two dif- 
ferent threshold hamming distances. For the smaller value of 
/ithr, V relaxes to a steady state value. At the larger value, 
the disease does not persist, (b) The hamming distance from 
the founder strain averaged over all strains present, (c) The 
standard deviation of the distance from the founder strain. 
This is a measure of the diversity of strains present. 



observe three long term behaviors. The first is trivial: 
when z < 1, the size of the epidemic goes to zero since on 
average each individual will expose fewer than one other 
individual. The remaining behaviors are (1) the fraction 
of the population that is infected, p, approaches either a 
steady state nonzero value or (2) after a brief outbreak, 
the epidemic dies out. Figure [l]a shows the difference in 
these two types of behavior. Interestingly, the transition 
from persistence to extinction occurs as we increase z. 

For a more complete characterization of this route to 
pathogen extinction, we ran simulations at every integer- 
valued point in the z-hthr plane with < /ithr < 10 and 
< z < 30. Figure 2 shows the result: below the curve, 
the epidemic persists, and above the epidemic dies out. 
Interestingly, /i t hr is a decreasing function of z, mean- 
ing a high contact rate is actually detrimental to the 
pathogen's ability to persist in a population. We can 
explain this result intuitively: the greater the value of 
z, the greater the average number of previous infections 
any given host will have. This implies that every individ- 
ual will have a larger memory repertoire, and thus each 
strain's ability to infect a host will be reduced. 

Increasing /i t hr also causes extinction to occur. This 
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FIG. 2. Boundary between pathogen persistence and ex- 
tinction in the bitstring model. Results were obtained by 
averaging the greatest value of h t hi at which the pathogen 
persisted over 150 runs, for each value of z. 



happens because the number of strains to which an in- 
dividual is immune increases with ft, t hr , and thus reduces 
the likelihood of infection. Increasing /i t hr can also be 
thought of as decreasing the pathogen mutation rate, 
since the memory immune response will be more effec- 
tive if the pathogen changes less rapidly. 

When persistence occurs, the fraction of the popula- 
tion infected is close to one (Figure pt). This is in con- 
trast to results from single strain SIR models, which typ- 
ically predict that the fraction of infected individuals is 
proportional to (and less than) 1 — 1/ Rq, where Rq is the 
number of new infections that would result from a single 
infection in a completely susceptible population [Q - this 
is z in our case. In Figure [lpt, the fraction infected at 
equilibrium is approximately 0.99, certainly greater than 
1 — 1/z = 0.8. This discrepancy occurs because there is 
more than one strain present in the bitstring model, and 
individuals can be infected by two different strains in two 
consecutive time steps. Figures |l]b-c illustrate this effect: 
in the case when persistence occurs, the mean hamming 
distance from the founder increases steadily in time, and 
the diversity of strains present (measured by the stan- 
dard deviation) is high. In contrast, when z is increased 
and extinction occurs, one can see that although the di- 
versity is initially higher, prior to extinction it is lower 
than in the persisting case. 

In the absence of immunity, no single strain has a com- 
petetive advantage over another in the bitstring model. 
However, there is evidence to suggest that in the case 
of influenza, there are selective pressures that constrain 
base pair subsititions to a small part of the entire genomic 
sequence space (i.e., with little diversity) and the num- 
ber of substitutions increases linearly in time || (Fig. |^). 
This motivates us to consider a special form of the model 
where bitstring mutation goes only in one direction, e.g., 
000. ..000 -> 100.. .000 -> 1100. ..000, etc. Specifically, 
when mutation occurs, rather than flipping a randomly 
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FIG. 3. Nucleotide substitutions in the nonstructural genes 
of influenza A. Figure reproduced with permision from Buon- 
agurio et al. f9[ 



chosen bit, the leftmost zero in the bitstring sequence 
is changed to a 1. All other model mechanisms are as 
before. 

Figure || shows that under these assumptions, the sys- 
tem's dynamics are quite different. Most notably, as the 
contact rate z increases from z — 1.25 to z = 2, sustained 
oscillations emerge from what appears to be steady state 
behavior surrounded with stochastic noise. As these os- 
cillations increase in amplitude, the minimum number 
of infected individuals ultimately gets so low that the 
epidemic dies out due to the population's finite size (at 
2 = 4). 

Figure [5| depicts the dynamics in sequence space for 
the last two cases from Fig. |j. In contrast to Figure [l], 
the hamming distance from the founder strain increases 
linearly in time and the diversity is low. These results 
mirror what has been observed for influenza, and mo- 
tivate us to understand the dynamics of this particular 
system in more detail. 

Under the assumption that pathogen evolution is con- 
strained to be linear in time, a further simplification of 
the system may be obtained by assuming that individuals 
are immune for a fixed period of time after infection. In 
the following sections, we analyze a difference equation 
model for this scenario to gain insight into the series of 
transitions observed in Figure |4| 
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FIG. 4. The dynamics of infection when bitstring evolution 
is constrained to go in one direction. We have used h t hr ~ 4. 
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III. IMMUNITY MODEL 

A. Model Derivation 

In this section, we derive a system of difference equa- 
tions to model the average behavior of a closed popu- 
lation of susceptible, infective, and immune individuals. 
Upon infection, individuals spend one time step in the 
infective class, and the subsequent r — 1 time steps in 
the immune class. After these r time steps, individuals 
are returned to the susceptible pool. The total popula- 
tion size is N and there are no births or deaths (i.e., N 
is constant). 

We define pt+i as the probability that an individual is 
infected at time t+1: 



Pt+i = x t s u 



(1) 



where Xt is the probability that an individual is exposed 
at time t, and St is the probability of residing in the 
susceptible class. At each time step, the Np t infected 
individuals make z random exposures. The probability 
that an individual is not involved in any such encounter is 
simply (1 — l/N) NptZ , and so the probability of exposure 
can be expressed as: 



^ \ Np t z 

n) 



(2) 



The fraction of the population that is immune at any 
given time i, can be determined by examining the fraction 
which has been infected in any of the previous r — 1 time 
steps, 521=1 Pt-k- The probability that an individual 



FIG. 5. The dynamics of infection in the bitstring model 
with directed evolution, (a) Fraction of population infected 
versus time for two different threshold hamming distances. 
For the smaller value of ftthr, p converges to oscillatory be- 
havior. At the larger value, the disease does not persist; the 
minimum number of infected individuals gets too low. (b) 
The hamming distance from the founder strain averaged over 
all strains present. The distance increases linearly with time, 
in contrast to the results in Figure |l|. (c) The standard devia- 
tion of the distance from the founder strain. This is a measure 
of the diversity of strains present. 



is susceptible is simply the probability of being neither 
infected nor immune: 



St = 



T-l 

1 - ^Pt-k- 

k=0 




Np t z^ 



T-l 



(3) 



(4) 



fe=0 



-e-«*)(l-5>_ fc ) (5) 

k=0 

in the large system limit, N — > oo. Thus, we have re- 
duced the immunity model to a r dimensional map, or 
equivalently a system of r difference equations. For large 
populations, the equations are independent of N, and 
thus the only parameters are r and z. This model was 
originally introduced by Cooke et al. M. 
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FIG. 6. Fraction of population infected versus time for r=6 
and various values of z. (a) For z = 2, p t relaxes to a fixed 
point, (b) For slightly larger values of z, pt exhibits small am- 
plitude quasiperiodic oscillations, (c) At z = 11, true period 
15 behavior emerges, (d) For z — 13, p t exhibits quasiperiodic 
oscillations that are heavily weighted by small values. 
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FIG. 7. Fraction of population infected versus time for 
r = 6 and z = 13, plotted on a semilog scale. The data 
set used in the above plot is the same as the one used in 
figure ^d. 
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B. Dynamics 

Numerical iteration of Eq. (j^) from the initial condi- 
tions pi = 0, i = l...r — 1 and p T = 10~ 4 yields the same 
four types of long term behavior observed in the bitstring 
model: approach to the trivial equilibrium (p t — > 0, not 
shown), approach to a nonzero equilibrium (Fig. |6|a), 
sustained oscillations that are generally quasiperiodic 
(Fig. ||b-c), and dynamic extinction (Fig. ||d). Dynamic 
extinction occurs only due to numerical roundoff error, 
and is not an analytical feature of Eq. (^) . As we increase 
z or r, pt eventually becomes so small at the minimum 
of the oscillation that (1 — e~ zpt ) in Eq. (||) numerically 
evaluates to zero. Furthermore, oscillations do not occur 
for r < 2. 

For a fixed r, we observe that pt relaxes to a fixed point 
for small enough z. When z < 1, p t approaches the trivial 
zero solution. As z is increased through z = 1, the non- 
zero equilibrium becomes an attractor. At some larger 
value of z, the fixed point loses stability and small am- 
plitude quasiperiodic oscillations appear symmetrically 
centered around the former equilibrium point. As z is 
further increased, the oscillations grow in amplitude and 
the system spends a large fraction of its time with only 
a small portion of individuals infected. Careful examina- 
tion of the oscillations in this regime suggest they consist 
of two phases: exponential growth followed by rapid de- 
cay (Fig. 

For certain combinations of r and z, the oscillations 



FIG. 8. The boundary between persistence and extinction 
for different population sizes. For each value of r, z was in- 
creased until the minimum expected fraction of infected indi- 
viduals declined below l/N, for N = 10 5 , 10 10 , and 10 15 . The 
epidemic is considered extinct in the upper right region of the 
figure because the minimum of the oscillation is below 1 /N . 



are truly periodic rather than quasiperiodic. We observe 
few patterns in the location of these periodic solutions 
in the parameter space. However, we note that in the 
somewhat rare cases in which periodic behavior emerges, 
the period is generally longer for larger values of z and 

T. 

In the following sections we use linear stability analy- 
sis to determine the location of the transitions from the 
trivial equilibrium to the nonzero steady state and from 
the nonzero steady state to oscillations. The transition to 
dynamic extinction is determined by deriving an approxi- 
mate expression for the minimum value of the oscillation 
and postulating that extinction occurs when this value 
goes below l/N where is any desired population size. 
Figure || illustrates the location of these transitions in 
the z — t plane. 



C. Stability Analysis 

The fixed points of the system satisfy 
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(l-e-»)(l-Tp*). 



(6) 



Making the substition q t — pt — p* and linearizing about 
P*, we obtain 



qt+i = (ze-^{l-Tp*) + e- z ^ -l) 



T-l 



-(e-^-ljE*-*- 



(7) 



fe=0 



Introducing the eigensolution q t — q \ l yields a polyno- 
mial for the roots A^: 

A r + «A r - 1 +/3A r - 2 + ...+/3A + /3 = 0, (8) 

where a = 1 — e~ zp (1 + z — zrp*), and (3 = 1 — e~ zp . 

When the t roots to this equation all lie within the unit 
circle in the complex plane, the fixed point in question 
will be stable. 

Case i. p* = 0. The first solution to Eq. (ra) is 



P 



0. 



At this point, the eigenvalue equation simplifies to 
A T - z\ T - 1 = 0. 



(9) 



(10) 



The first r — 1 roots of this equation are A = 0. The final 
solution solves A — z — 0, and thus the only eigenvalue of 
interest for stability is 



A = 



(11) 



indicating that the trivial equilbrium is stable for all 
z < 1. Not surprisingly, this agrees with our previous 
results from the bitstring model in Section [Hj, since if 
z < 1, fewer than one new infection will result from each 
currently infected individual. 

Case ii. p* ^ 0. The results in Figure |] suggest that 
when z > 1, the nonzero solution to Eq. (^) becomes 
stable. Indeed, Cooke et at have shown for all z > 1 
that this point is globally stable when t = 1 , and locally 
stable when r — 2. They conjectured that when r = 3, 
the fixed point loses stability at z = 4.58. Here we will 
verify this result and obtain the transition for all t > 3 
by finding the location of the bifurcation at which the 
quasiperiodic orbits emerge. 

As before, the onset of instability occurs when |A»| = 1 
for one or more A^. Therefore, to locate the transition, 
we substitute A = e 4 * into Eq. (|J). This yields two new 
equations, one for the real part and one for the imaginary 
part: 



1 



-zp 



— (1 — Tp*)ze~ 



zp _ 



cos(0) — cos(0(r + 1)) 



1 



cos(t</>) — 1 
2 sin(</>) - sin(2</>) 



sin(</>(r — 1)) + sin(0) — sin(r</ 



(12) 



(13) 



Combining these two equations with the expression for 
p* (Eq. (ph), we have three equations for four unknowns 
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FIG. 9. Hopf curve generated by numerically solving equa- 
tions (^|) and ( [I4 ) simultaneously for specified values of r. The 
open circle marks the location of the transition previously de- 
rived by Cooke et al. M. The solid line is our asymptotic 
solution, taken to 0(1/-P). 



(p*,T,z, and <f)). These three equations define the Hopf 
curve. 

Figure || shows the numerically generated solution for 
the Hopf curve in the z-t plane. The dependence on 
parameters is clearly similar to what we observed in the 
biststring model as that system switched from persistence 
to extinction. For z < 1, p* = is the only attractor. 
When z > 1 below the Hopf curve, the stable fixed point 
is given by the non-zero solution to Eq. (^). 

For large r, we can write a perturbative solution for 
the Hopf curve. Defining a new variable e = we can 
express <f>, z, and p* as power series in e: 



z = 1 



00 



(14) 
(15) 
(16) 



Solving perturbatively for the ai, bi, and Cj gives 
* = 1 + y e + ^ + A ( i2^ + ,4 )£ 3 + 0{e% (17) 



= 7r e+|e 2 + Je 3 + O(e 4 ). 



(18) 
(19) 



Figure ^ shows that for large r, the perturbative solution 
for the bifurcation curve agrees well with the numerical 
solution. 

The variable <fi may be thought of as the rotation num- 
ber |l(J of the solution to the linearized equations. In 
cases where a periodic orbit emerges at the bifurcation, 
27r/0 will be the period of the orbit, and in cases where 
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FIG. 10. Generic form of a relaxation oscillation for large z 
and r. The oscillation is characterized by exponential growth 
followed by rapid decay. 



the orbit is quasiperiodic, the approximate pattern will 
repeat itself on average every 2n/(p time steps. 

Eq. ( |l9| ) indicates that tf> is a decreasing function of 
t = 1/e, meaning that the period of the epidemic oscil- 
lations will increase with r. In other words, we expect 
oscillations to occur less frequently as the duration of 
immunity increases. 

D. Relaxation Oscillations and the Route to 
Extinction 

As pointed out earlier, Figure suggests that as z and 
r become large, the system's dynamics can be character- 
ized by exponential growth followed by rapid decay, with 
short transitional phases between the two regimes. In 
this section, we explore the relaxation oscillations by de- 
riving equations that approximate the system's behavior 
in the various phases. 

We begin by focusing on a single oscillation as pic- 
tured in Figure [l^. As pt grows exponentially toward 
its maximum, the fraction of individuals infected at the 
current time and the previous r — 1 time steps is small 
enough that, to a first order approximation, Eq. (^) can 
be rewritten as 

Pt+i ~ zp t (in the growth phase). (20) 

The behavior persists until zpt reaches order 1. We make 
the assumption that Eq. ( po| ) holds until the point po and 
check for consistency after subsequent calculations. For 
times less than t — 0, we can write p-t = Po/ z f . 

Next, we iterate the r-dimensional map using the un- 
kown form of po to find pmin{po)- The tilde notation in- 
dicates that p m in is a local minimum of the oscillations. 
To find the global minimum, we must find the value of 
Po that minimizes p m i n . In what follows, we derive an 
approximation for p m i n in the large z limit. This calcu- 
lation requires that the points near the transition phase 
be handled individually before a general expression for 
the decay behavior may be obtained. We determine pi, 
P2, and P3 explicitly and then derive the general form for 
Pt+i for t > 3. 



A reasonable approximation for p\ — xqSq is obtained 
by including only the first term that appears in the sum 
in s : 

Pi«(l-exp(-zpo))(l-po) (21) 

Next, p2 can be determined by first considering the 
probability of exposure at time t = 1. Since the fraction 
of individuals infected is order one at the maximum point 
of the oscillation, the probability of exposure at time t = 
1 goes to one in the limit of large z. This implies that 
virtually all individuals who are susceptible at time t = 1 
will become infected at time t — 2: p2 ~ s\. 

In order to produce a simple expression for p2, it useful 
to note that the fraction of individuals which reside in the 
susceptible class at any point in time can be written in 
terms of the same fraction at the previous time step: 

s t = exp(-zp t -i)s t -i + pt-r- (22) 
Using this, we write P2 in a convenient form: 

p 2 « si w exp(-zp )(l -Po) + (23) 

To find P3 , the final point in the transition region, we 
note that S2 ~ P2-t ~ since zp\ 3> 1. Assuming 

zp2 is small, we obtain: 

P3 ~ zp 2 s 2 = ~^P2- (24) 

The general behavior in the decay regime can be de- 
rived by examining the fraction of individuals susceptible 
for 3 < t < t. In this region, St-i is small compared to 
Pt-r, yielding: 

at « Pt-r « -^=- t (for 3 < t < r). (25) 

Since zp t is small in the decay phase, the probability of 
exposure depends linearly on the fraction of individuals 
infected: x% « zp t . Thus, we see that for t > 3, we can 
again replace the original set of r difference equations by 
a single equation: 

Pt+i = z'^PoPt- (for 3 < t < t). (26) 

The minimum of p t occurs when z t ~ T+1 po becomes 
greater than 1. If we assume that po lies between - and 1, 
as is consistent with our earlier assumption that Eq. (|2p| ) 
holds until t = 0, then the minimum must occur at t = r. 
Consequently, we can express pt in terms of p3 for times 
greater than t — 3. Combining Eqs. ( p6| ) and (p3|), we 
obtain: 

Pr=z-^ T2 - 5T+6) P T - 2 P2 (27) 

Finding the minimum value of p T requires solving for 
the roots of the equation dp T /dpo — 0, which yields 
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+P | — — 1 -Po) J =0 



We try a solution of the form 

log(/) 
Pa = • 



(28) 



(29) 



Substituting this into Eq. (|28|) we obtain an asymptotic 
expression for /, which is valid to leading order as z — > 
oo: / « z T /(r — 1). This gives: 



Po 



rlog(z) 



(30) 



which is consistent with our assumption that po must lie 
between 1/z and 1. 

Finally, to obtain p m i n , we insert (|30[ ) into our expres- 
sions for p 2 and p T (Eqs. (|23|) and (|27|)): 



„ . - z -|(r 2 -r+2) „^T-1 
mm — * 



(rlog(z)r- 1 . (31) 
In the derivation of the difference equations in Section 



III A , we assumed an infinite population size. Under this 
assumption, the fraction of immune individuals can be 
arbitrarily close to 1 without driving the disease to ex- 
tinction. If, however, we assume the population size is 
finite, at some point the minimum value of pt will be less 
than the fraction of the population equivalent to one in- 
dividual, that isp m i„ < 1/N. Replacing p min by 1/N, we 
have an equation for the curve that separates the region 
of disease persistence from the region of extinction: 



1 

N 



r+2) 



(rlog(z)) T ~ 



(32) 



Figure (11) compares the predictions of ( |32| ) with the 
data obtained from the map dynamics. We see that our 
predictions work quite well for large z. Note that the 
asymptotic approximation only fits the numerical data 
for unrealistically large population sizes (for smaller pop- 
ulation sizes the extinction boundary occurs at smaller 
z). This is not too troublesome, however, since our 
overly-simplified model cannot be expected to quanti- 
tatively match actual population features. Rather, the 
strength of this approach is that it provides a clear pic- 
ture of a mechanism for dynamic extinction. 



IV. CONCLUSIONS 

We have shown that oscillations in the number of in- 
fected individuals in a population could be due to a mu- 
tating pathogen. In both of the models we have stud- 
ied, oscillations occur as a consequence of the continual 
introduction of novel strains, rather than the interplay 
between several pre-existing strain types (for studies of 
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FIG. 11. Extinction curves for various values of system size 
TV. The solid lines represent the asymptotic approximation 
for large z. Here, TV is merely a parameter that determines 
the condition of extinction for the difference equation model; 
it is not the number of nodes in a simulation as in the bit- 
string model. The disease is considered extinct if the fraction 
infected ever falls below 1/N. 



the latter, see references [TlJ and p2j ) . We can explain 
this phenomenon naturally in terms of single outbreak 
epidemics: each period of oscillation can be regarded as 
a new epidemic with a new strain to which few if any 
individuals have immunity. 

In the r-dimensional map, oscillations occur only for 
an intermediate range of immunity duration r. If it is too 
short (or, equivalcntly, if the mutation rate is too high), 
oscillations do not occur because a significant pool of 
susceptible individuals is always present. Alternatively, if 
the duration of immunity is too long, the infected pool os- 
cillates with increasingly large amplitude and ultimately 
becomes too small at its minimum for the pathogen to 
persist. 

The presence of oscillations depends on the contact 
rate in a similar fashion. For low contact rates, the sus- 
ceptible pool remains large and oscillations do not occur. 
For high contact rates, large amplitude oscillations force 
the number of infected individuals to such a small value 
that the epidemic dies out. 

We have quantified the effect of contact rate and im- 
munity duration on oscillatory behavior through Hopf 
bifurcation analysis to locate the onset of oscillations, 
and with asymptotic methods to determine their mini- 
mum value. In both analyses, the boundary between two 
types of behavior is marked by an inverse relationship 
between immunity duration and contact rate. In partic- 
ular, in a population with high contact rates, pathogen 
persistence requires short periods of immunity (or high 
mutation rates) suggesting that highly connected popula- 
tion structure could provide selective pressure for rapidly 
mutating pathogens. This could lead to diseases that are 
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more difficult to actively counteract through immuniza- 
tion programs. 

There are two assumptions that could have significant 
impact on our results. First, these models operate in dis- 
crete time, where oscillations and chaotic dynamics are 
known to occur more readily. Second, we have assumed 
the underlying networks are fully mixed, but some sort of 
quasi-static network structure may be much more realis- 
tic. The addition of spatial effects could dampen oscilla- 
tions, if two regions oscillated out of phase and thereby 
prevent global extinction of the pathogen. 
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